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(Nj , Abstract 

£ : 

c3 . We propose a locally one dimensional (LOD) finite difference method for multidimensional 
^~ "i • 

■ Riesz fractional diffusion equation with variable coefficients on a finite domain. The 
numerical method is second-order convergent in both space and time directions, and its 
unconditional stability is strictly proved. Comparing with the popular first-order finite 
difference method for fractional operator, the form of obtained matrix algebraic equation 
is changed from (I - A)u k+1 = u k + b k+1 to (I - A)u k+1 = (I + B)u k + ~b k+1 ' 2 - the 
three matrices A, A and B are all Toeplitz-like, i.e., they have completely same structure 
and the computational count for matrix vector multiplication is O(NlogN); and the 

computational costs for solving the two matrix algebraic equations are almost the same. 

> ■ 

(■T) ■ The LOD-multigrid method is used to solve the resulting matrix algebraic equation, and 
\q ■ the computational count is O(NlogN) and the required storage is O(N), where N is the 



number of grid points. Finally, the extensive numerical experiments are performed to 



show the powerfulness of the second-order scheme and the LOD-multigrid method. 
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1. Introduction 

In recent years considerable interests in fractional calculus have been stimulated by 
the applications in physical, chemical, biological, and engineering, etc., areas jsj]. The def- 
initions of fractional calculus are versatile, e.g., Riemann-Liouville derivative, Grunwalc - 
Letnikov derivative, Caputo derivative, Weyl derivative and Riesz derivative et al fill. [l4|. 
and they are not completely equivalent. Depending on the particular applied field, some- 
times one of its definitions is more popular than others. For example, the Riesz fractional 
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derivative appears in the continuous limit of lattice models with long-range interactions 



19j. This paper focuses on the multidimensional Riesz fractional diffusion equations. 

Nowadays, the finite difference discretization for space fractional derivatives is expe- 
riencing rapid development, including the Riesz fractional derivative; such as, Yang et 
al numerically study the Riesz space fractional PDEs with two different fractional orders 
1 < a < 2 and < (3 < 1 24[; Zhuang et al consider a variable-order fractional advection- 
diffusion equation with a nonlinear source term on a finite domain 26[. In the last two 
years, for the space fractional derivatives, we notice that two different second-order dis- 
cretization schemes are developed {17I . [20 ] ; even the third-order discrezation scheme is 
obtained |25l| if a compact difference operator is performed on the discrezation scheme 
given in [20|. 

Another topic related to effectively solving the equations involving fractional operators 
is about how to efficiently solve the resulting matrix algebraic equations. The 'unlucky' 
thing is that the matrix in the matrix algebraic equation is usually full because of the 



nonloca 



12 
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properties of the fractional operators, and the 'lucky' thing, as pointed out in 



22[, is that the matrix has some special structure, i.e., the matrix is Toeplitz- 



like matrix, and the count of its matrix vector multiplication is O(NlogN) by using 
the constructed circulant matrix and fast Fourier transform, and the required storage 
is O(N), where N is the number of grid points. Pang and Sun 12[ successfully use 
the multigrid method (MGM) to efficiently solve the resulting matrix algebraic equation 
of the one dimensional fractional diffusion equation by using first-order discretization 
scheme 10j. Here we further extend the MGM to solve the matrix algebraic equation of 



the multidimensional Riesz fractional diffusion equation discretized by the second-order 
scheme. 

We use the LOD strategy to solve the multidimensional Riesz fractional diffusion 
equation. The LOD methods include alternating direction (AD) methods and fractional 
step procedures [6|. The AD methods were first introduced in three papers 0, Q, H] 
by Douglas, Peaceman, and Rachford. The Peaceman and Rachford (PR-AD) method 
works well for two-dimensional problems. However, it can not be extended to higher 
dimensional problems. Douglas (D-AD) method are valid for any dimensional 

equations. And PR-AD and D-AD are equivalent in two-dimensional problems. Both PR- 
AD and D-AD schemes are used in this paper to discretize the multidimensional Riesz 
fractional diffusion equation. In each dimension, the obtained matrix algebraic equation is 
solved by MGM. Although the spacial fractional derivative is discretized by second-order 
scheme, for any single dimension the form of the obtained matrix algebraic equation is 
(I — A)u k+1 = (I + B)u k + £r +1 / 2 , in fact the corresponding form for the first-order 
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discretization scheme is (/ — A)u k+1 = u k + b k+1 ; the three matrices A, A and B are all 
Toeplitz-like, and they have completely same structure and the computational count for 
matrix vector multiplication is O(NlogN); and the computational costs for solving the 
two matrix algebraic equations are almost the same. In other words, the second-order 
scheme improves the accuracy but almost without increasing the computational cost. 

More concretely, in this paper using the second-order accurate and unconditionally 
stable computational scheme and LOD-MGM, we solve the following multidimensional 
variable coefficients Riesz fractional diffusion equation with the computational count 
O(MogiV) and storage 0(N), 

' du(x,y,z,t) d a u(x,y,z,t) d^u(x,y, z,t) 

— m — = c(X; v " z > t] Wr (x ' y ' z ' } dw 

d^u{x,y,z,t) 
+ e(x,y,z,t) + f[x,y,z,t), 

u(x,y,z,0) = u (x,y,z) for (x,y,z)e£l, 
u(x,y,z,t) — for (x, y, z, t) £ dVL x (0, T], 

in the domain fl = (xl, xr) x (yz, ya) x (zl, zr), < t < T, where the orders of the Riesz 
fractional derivatives are 1 < a,/3,7 < 2; f(x,y,z,t) is a source term and the variable 
coefficients c(x, y, z, t) > 0, d(x, y, z, t) > O^eix, y, z, t) > 0; the Riesz fractional derivative 
for n E N, n — 1 < v < n, is defined as"0, Esl 

d u u(x, y, z, t) / ,, \ , x , 
^ ; = -k v { XL D»+*D: R )u(x,y,z,t), (1.2) 

where the coefficient k v = t, — r — t^, and 

" 2cos(v-k /2) ' 

x L D v x u(x, y, z, t) = — - ^ — J (x-CT'^ui^y^,^, (1.3) 
(—]) n Ff 1 r x R 

x D» XR u(x, y, z, t) = Jl-J—— J (£ - x)"-"-^^ I/, ^ (1-4) 

are the left and right Riemann-Liouville space fractional derivatives, respectively. 

The outline of this paper is as follows. In the next section, we introduce the second- 
order finite difference discretizations for the Riesz fractional derivatives; and the full 
discretization of fll.lD is derived, where the Crank-Nicolson scheme and LOD method 
are combined together. We theoretically prove the presented finite difference scheme is 
unconditionally stable in Section 3. In Section 4 we propose a V-cycle LOD-MGM for 
the resulting system of (II. ip . To show the powerfulness of the second-order scheme and 
LOD-MGM, the extensive numerical experiments are performed in Section 5. Finally, we 
conclude the paper with some remarks in the last section. 
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2. Derivation of the finite difference scheme 



In this section, we derive the full discretization schemes of fll.ip . The first subsection 
introduces the second-order finite difference discretizations for the Riesz fractional deriva- 
tives in a finite domain. Then in the second subsection, we present the scheme for the 
one dimensional case of fll.ip . The third and fourth subsections detailedly provide the 
two dimensional case of ( II. ip and (11.11) itself, respectively. 

2.1. Discretizations for the Riesz fractional derivatives 

Take the mesh points Xj = Xl + iAx, i — 0,1, ... , N x , i/j — y L + jAy, j — 0,1, ... , N y , 
z\ = zl + IAz, I = 0, 1, . . . , N z and t k = kAt, k = 0, 1, . . . , N t , where Ax = (xr — Xl)/N x , 
Ay = {yR — yhj/Ny, Az = (zr — zl)/N z , At = T/N t , i.e., Ax, Ay and Az are the uniform 
space stepsizes in the corresponding directions, At the time stepsize. For v G (1,2), the 
left and right Riemann-Liouville space fractional derivatives ( II. 3p and (II. 4p have the 
second-order approximation operators S^ +X u^j tl and <J V) _ z it* . j, respectively, given in a 
finite domain (3,[l7|, where it* - j denotes the approximated value of u(xi,yj, zi,t k ). 

The approximation operator of (II. 3p is defined by 0, 17] 



t+i 



Si/, + X U i j l '•' 



T(4- v){AxY 



Ev k 



m=0 



and there exists 



where 



XL D v x u(x, y, z, t) = 5„ i+x ul jtl + 0(Axf 



(2.1; 



(2.2) 



1, m = 0, 

-4 + 2 3 -», m = 1, 

6-2 5 ^ + 3 3 ^, m = 2, 
(m + lf~ u - 4m 3 - u + 6(m - l) 3 ^ 

- 4(m - 2f~ u + (m - ?>f~ y , m > 3. 

Analogously, the approximation operator of fll.4j) is described as 

, N x -i+l 



(2.3) 



T(4 - i/) (Ax 



- V o v v k 



m=Q 



and it holds that 



,D v XR u{x, y, z, t) = 5 Vt _ x u^ + 0(A 



x 



(2.4) 



(2.5) 



where g v m is defined by ( 12. 3ft . 

Combining (12.21) and (12.51) . we obtain the approximation operator of the (Riemann- 
Liouville) Riesz fractional derivative 
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d u u(x,yj, z h t k ) 



d\x\ u 



K v (xiPZ +x D x R ) u(x, y s , z h t h 
k v (6 Vi+x + S Vi _ x ) ul jtl + O(Ax) 2 

t+l N x -i+l 

E v f. 
9rn U i+m-l,j 



T(A-u)Ax» (E 5 



m i— m+1,3 



-m+l u m 



r(4 - u)Ax 



f- 0(Ax) 2 

,v + E ^-m<^ ) + ( Aa; ) 2 



m=0 

AT* 



^m=0 
iV. r 



m=i—l 
2 



(2.6) 



where 



dim 



9i—m+\i 


m < % — 1, 


g u + 92, 


m — i — 1, 




m = i, 


9% + 9%, 


m — i + 1, 


k flm— i+lJ 


m > i + 1, 



(2.7) 



with z = 1, . . . , iVa — 1, together with the Dirichlet boundary conditions that define «ojz 
and u k Nx j^ as appropriate. 

Taking z/ = 2, both Eq. (12.21) and (12. 5p reduce to the following form 

d 2 u(xj,y,z,t) _ u(x i+1 ,y, z, t) - 2u(x h y, z, t) + ttfa-i, y, z, t) 2 

dx 2 ~ (Ax) 2 + 1 } ■ 

Similarly, it is easy to get the one- dimensional and two-dimensioanl case of ( 12.ip -( l2T7j) . 

2.2. Numerical scheme for ID 

Consider the one- dimensional Riesz fractional diffusion equation 



du(x,t) , . <9 a w(x,t) . 



Of 



d\x\ 



(2- 



In the time direction, we use the Crank-Nicolson scheme. Taking the uniform time 
step At and space step Ax, and setting c\ = c(x i: tk) and f^ +1 ^ 2 = /(xj, ^+1/2), where 
tk+1/2 = (tk + tfc + i)/2, the full discretization of (12. 8 p has the following form 



u h+1 - u k 



At T(4 - a) Ax 



fc+1/2 n x k , k+l 

\ " ~™ «m "I" U m rk+1/2 



/ j Hi,n 
m=0 



(2.9) 
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Then ( 12. 9 j) can be expressed as 

At 



1--51 U^=[l + ^CK + A ^ , 



c fc+l/2 



(2.10) 



where 



_ fc+l/2 N x 



T(4 - ct)Ax a 



Putting ^f +1//2 = 2 r(A-a)Ax c ' ' *" ne s y s ^ em °f equations given by (I2.10p takes the form 



m=0 

fc+l/2 



T(4-ct)Ax 



fc+l/2 % 



m=0 



(/ - A k+1/2 )U k+1 = (/ + A fc+1 / 2 )f/ fc + AtF fc+1 



/2 



(2.11) 



where I is the (N x — 1) x (N x — 1) identity matrix, 

U k = [u$, «*..., <„J T , F fc +V 2 = ^ *+V> / . +1 / 2) _ _ _ ^ / ,+i/2 ]T) 

and the discretizations at the interior x-gridpoints define the entries of the matrix A k+1 ' 2 , 
A k ^ 2 for % = 1, . . . , N x — 1 and m = 1, . . . , N x — 1 are defined by 



.4 



fc+l/2 



fc+l/2 



a x c fc+l/2 



a t" 
Hi— m+lS i 

+ ^)ef +1/2 , 

a .fc+l/2 



m < z — 1, 
m — i — 1, 
m — i, 
m = z + 1, 
m > z + 1. 



(2.12) 



2.5. I CD scheme for 2D 

Consider the following two-dimensional Riesz fractional diffusion equation 

^ ^ - c(x, y, t) + .(x, y , t) ^fe^ + /( x, y , t). 



Ot 



d\x\ 



d\y 



(2.13) 



Analogously we still use the Crank-Nicolson scheme to do the discretization in time 
direction. Taking • as the approximated value of u(xi,yj,tk), c k j = c(xi,yj,t k ), d k j = 
d(x i7 yj,t k ), t n+1/2 = (t„ + t n+l )/2, ftj~ 1/2 = f(xi,yj,t k+1 / 2 ), Ax = (x R - x L )/N x , and 
Ay = (yR — yi)/Ny, for the uniform space steps Ax, Ay and the time stepsize At, the 
resulting discretization of (12. 13[) can be written as 



„ k + l _ _k r k+l/2 Nx v k + v k+1 

hi hi a hi \ ^ —a m,i ' m,j 



At T(4 - a)Ax Q ■ 

—Kpd\ 



^1 9i,r 



+ 



m=0 

fc+l/2 Ny 

%3 



r(4-/3)Ax/ 3 ^^' m 

v ' ' m=0 



7 ,fc _i_ 

E—/3 "i.m ^ "i,m rk+1/2 
9j,m rt ' J i. 



(2.14) 



Similarly, we define 

fc+1/2 TV, fc+1/2 ATj, 

5" „*. = ^iiJ v„« u fc .. s" u k+1 = V? u k+1 - 

^ ' m=0 ^ ' m=0 



X» „fc V" ~P k en fc+1 _ V ~/3 fc+1. 

U _ (3) Ay 13 Z—i y J' m *' m ' P>v tj y{A - (3) Ay 13 ^ j ' m ,m ' 

then Eq. (I2.14p can be rewritten as 

i - ft - ft) «S* = (i + t C + + At ^" 2 ' 

For the two-dimensional Riesz fractional diffusion equation (12.13jl . the relevant per- 
turbation of (12. f 5p is of the form 

yc) (i - f*,) «?.r = (i + f c) (i + t^) <, + 

(2.16) 

The scheme (12. 16[) differs from (12 . 1 5[) by a perturbation 



^—Ls" X" (v k+1 - v k \ 

which may be deduced by distributing the operator products in (I2.16p . Since {u^ 1 —u k ^) is 
an O(At) term, it follows that the perturbation contributes an 0((At) 2 ) error component 
to the truncation error of (12.151) . Thus, the scheme (12.16)) has a truncation error also 
0{{AxY) + 0{{Ayf) + 0{{Atf). 

For efficiently solving system (12. 16 j) . the following techniques can be used: 

D-AD scheme 

1 " %%,v) <f = < - ^Af> (2-18) 



where w* ■ is an intermediate solution. Subtracting ( 12 . 1 8|) from ( 12 . 1 7j) . we obtain 



PR-AD scheme 



u i,j 



fc+ i At „ , k fc+1 x 



(l " f C) = (l + Y S 'L) <, + f (2 ' 19) 

(i - ' = f 1 + T «*) <' + T^" 2; (2 ' 20) 

with intermediate solution u* -. Subtracting (12.201) from (12.191) . we have 
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(2.21) 



2.4- D-AD scheme for 3D 

Similarly, the resulting discretization of ( II. ip can be written as, 

I 1 2 a ' x 2 fi,v 2 hZ ) iJ ' L 

fi^toj-,, At At \ fc .k+i/2 A . 
= ( 1 + y^v + -y ^ )2/ + -y V ) + h,j,i At. 

The perturbation equation of (I2.2ip is of the form 
The scheme (I2.22p differs from ( I2.2ip by the perturbation term 

(At) 2 , „ „ „ „ „ „ W fe+l_ fe N (At) 3 „ „ „ ( k+1 _ k \ 

^ \ a,x°/3,y + °q,x°7,2 + l3,y°-y,z)\ U i,j,l a i,j,U g °«,3; P,y y,z\ U i,j,l a i,j,l)- 

The system of the equations defined by (12.221) can be solved by the D-AD scheme (5, 6] 
(l - <i - (l + + + A*?,) <, + A^ (2-23) 



(2.22) 



(1 " ^) M S = <i " T^-'- (2 ' 25) 

For maintaining the consistency, we need to carefully specify the boundary conditions 
of u>™'ji and According to fl2.23p - f)2.25p . we obtain 

3. Convergence and Stability Analysis 

We show the convergence for one- dimensional and multidimensional Riesz fractional 
diffusion equation by proving the consistency and stability (according to Lax's equivalence 
theorem) . 

Lemma 3.1 (^). The coefficients g^ m , v G (1,2) defined in \2. 7| ) satisfy 



(1) ^<0, g» >m >0(m^i); 

N x N x 

( 2 ) Yl Km < and - > ^ g v i m . 

m=0 m=0,mj^i 
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3.1. The stability of the numerical methods in ID 

Theorem 3.2. The Crank-Nicholson scheme (12.111) of the Riesz fractional diffusion 
equation (12. 9p with 1 < a < 2 is unconditionally stable. 

Proof. First, we prove that the eigenvalues of the matrix A k+1 ^ 2 have negative real parts. 
Note that A^ 1 ^ 2 = <7fj£f +1//2 , and from Lemma 3.1 we obtain 

N x N x 
Ti= £ = -< 1/2 - (3-1) 

According to the Greschgorin theorem js], the eigenvalues of the matrix A k+l l 2 are in the 
disks centered at A k ^ 1 ^ 2 , with radius r i; i.e., the eigenvalues A of the matrix A k+1 l 2 satisfy 

\\-A k f 2 \<r h (3.2) 

thus, the eigenvalues of the matrix A k+1 / 2 have negative real parts. Similarly, we can 
prove that the eigenvalues of the matrix / — A k+l l 2 have a magnitude greater than 1 and 
invertible. 

Note that A is an eigenvalue of the matrix A k+1 l 2 if and only if 1 — A is an eigenvalue 
of the matrix / - A k+1 / 2 , if and only if (1 - A)-*(l + A) is an eigenvalue of the matrix 
(I-A k+l ' 2 y\I + A k+l l 2 ). Since 9fc(A) < 0, it implies that |(1 - A)~ x (l + A)| < 1. Hence, 
the spectral radius of the matrix (J — A k+1 ' 2 )~ l {I + A k+1 ^ 2 ) is less than 1. □ 

3.2. The stability of the numerical methods in 2D 

Under a commutativity assumption for the operators (l — ^S^x) and (l — ^#g i2/ ) in 
(I2.15p . the PR-AD scheme and D-AD scheme will be shown to be unconditionally stable. 
The commutativity assumption for these two operators is a common practice in estab- 
lishing stability of the classical AD methods for the diffusion • The commutativity 
of these operators implies that the matrices Ajdai anc ^ ^w^y gi yen m A3.7[) commute. 
Theorem 3.3. Both the D-AD scheme (l2T7jl - (12151) and PR-AD scheme (F2~T9l) - (l2~20j) . 
defined by (12 . 1 6j) . are unconditionally stable for a, (3 G (1,2), if the matrices A^ai an d 
^2£>Ar, commute. 

Proof. D-AD scheme (I2.17p - (l2.18p can be expressed in the form 

(I - 4 + n%)U* = (/ + J&% + 2A k 2 + D f y )U k + AtF k ^ 2 ; (3.3) 

(i - A + D %) uk+1 = u *~ 4if y u k ; (3.4) 

and PR- AD scheme (l2TT9|) -( l2T20|) is of the form 

(I - AttW = (/ + 4t%)U k + ^F fc+1/2 ; (3.5) 



9 



(I - 4 + nt) Uk+1 = ( J + A + d%) U * + ^F fc+1/2 ; (3.6) 
where the matrices -A^AAx an d Ar^Ay denote the operators ^<^£ an d ^&'p,y> an d 

f/ fc = [li^i, W 2 ,l) • • • 5 W^b-IjI) M l,25 u 2,25 • • • 5 ^A^-1,25 • • • 5 ^1,^-1) U 2,JV S -1; • • • 5 u ^ , ir -l l J^ / -l] T ) 

rrt r t * * * * * * # "iT 

U ~~ L"l,l> "2,1) • ' • 5 "JVas-1,1) "1,25 "2,25 • • • 5 "JV^-1,25 • • • 5 U l,Ny-\l u 2,N y -li • • • 5 " J\T a - 1 , JV V - 1 J 5 



and the vector F k+1 / 2 absorbs the source terms f i ■ and the Dirichlet boundary con- 
ditions at time t = tk+\ in the discretized equation. The matrices Ajo^ai an d A^J^ are 
matrices of size (N x - 1)(JV V - 1) X (N x - l){N y - 1). 

Let us cancel the intermediate solution U*, then D-AD scheme and PR-AD scheme 
have the same form 

(I ~ A^aIW - <Z> k+1 = (' + AnilW + 4Z> k + A * Fk+1/2 > (3-7) 
from (13.71) we have the perturbation equation 

(T A k+1 / 2 \(T _ (j , A k + 1 / 2 \(T , A k+l/2^r?k 

K 1 ~ A 2D,Ax)" ~ A 2D,Ay> h ~ K 1+ A 2D,Ax)\ I + A 2DAy> h ' 

where 

j?k _ r k k k k k k k k k iT 

^ — [^1,1 5 C 2,l5 • • • 5 C N x -l,ll C l,25 C 2,25 • • • 5 ^^-1,2) ' ' ' > e l,JV„-l5 e 2,JV„-l5 • • • 5 fi N x -l,Ny-l\ 5 

and e k j = u(xi,yj,tk) — u k j, consequently 

E k = [(/ - A\i%Y\l - A k +%y\l + A + D %)(I + A%fXE°. 
Since the matrices i^Ai an d ^dKj/ commute, it can be written as 

E k = [(/ - A k 2 + D %y\i + - a^)-^/ + 4 + D fX E °- 

According to Theorem 3.2, similarly, it is easy to check that the eigenvalues of the matrix 
^2DAx an d ^'dAj have negative real parts. Then, both the spectral radius of the matrixes 
( J - AId/ax)' 1 ^ + Ad Ax) and ( J - AdAv)' 1 ^ + Ad Ay) are less than X ' tnerefore the 



sequence 



BlfJ-^&r^ + ^i)]* a nd [(/-^^^(I + ^S^)]* converge to zero 



matrix 



16j. Hence, the difference scheme (I2.16P is unconditionally stable. □ 



3. 3. The stability of the numerical methods in 3D 

Theorem 3.4. The D-AD scheme (l2~23|) -( l2T25l ). defined by (|2~22l) . is unconditionally 
stable for a,/3, 7 G (1,2), if the matrices A 1 ^ 1 ^., A^ 1 ^ and A k ^^ z commute. 
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Proof. D-AD scheme fl2.23p - fl2.25p can be written as 



rrfc,l Ak+l/2 TT k. 
U ~ A 3D,Ay U ' 



fc+1/2 



, fc+1/2 
3D, Ay 



2A k 3 i 1 i 2 z )U k + AtF k+1 / 2 ] 



|Hl/2u r t,2 
l 3D,Ay) U 



according to (I3.8I) - (I3.10I) . we have the following equation 



A k+1/2 )(T 



A k+1/2 )(T 

/1 3D,AyJ l + 



fc+1/2 wi-fc+l 



A 2,D,Az) U 



/2 



(3.8) 
(3.9) 
(3.10) 

(3.11) 



where the matrices A^~£ x and A k ^^ y and A 1 ^ 1 ^ denote the operators ^S" x and 



2 "fty 



and 4^5" , respectively, the vector t/^' 1 and ?7 fc ' 2 denote the intermediate solution, the 
vector F k+1 / 2 absorbs the source terms /f^ 2 and the Dirichlet boundary conditions at 
time t = tfc+i in the discretized equation. The matrices A 1 ^ 1 ^ and A^ 1 ^ and A^ 1 ^ 
are matrices of size (N x - l)(N y - l)(N z - 1) x (iV a - l)(N y - 1)(N Z - 1). 

By the similar analysis, it can be proven that the spectral radius of the matrices 



(I - 4 + d%)-\I + (I - Al + D %r\I + AH%) and (/ - AH%)-\I ^ 

are less than 1, therefore the difference scheme (12.221) is unconditionally stable. 



fe+l/2 - 



fc+1/2 - 



fc+1/2 \ 
3D,Az) 

□ 



4. Multigrid method for the resulting matrix algebraic equations 

We use a V-cycle LOD-MGM to solve the resulting matrix algebraic equations of (11. ip . 
Meanwhile, we show the convergence of the resulting system. In order to develop a fast 
algorithm, i.e., realizing the computational count O(NlogN) and the required storage 



O(N), as did in 12|, |2l|, |22J, we first introduce the Toeplitz matrix and the circulant 
matrix. The n x n Toeplitz matrix T n (c) is defined by [l| 



T n (c) 



n 

j-fcjj,fc=l 



CO 



C-l 
CO 



C-(n-l) 
C-( n -2) 



(4.1) 



c n-l C n _2 ■ ■ • Co 

and the circulant matrices are the "periodic counsins" of Toeplitz matrices. We denote 
by circ (c , c±, . . . , c n _i) the circulant matrix whose first column is c = (c , ci, . . . , c n _i) T , 



Co 

Cl 
C'2 



Cn-l C n _2 
Co C n _i 
Cl Cq 



c 2 
c 3 



Cl 

c 2 
c 3 



Cn-2 C n _3 ' • " • Co C n _i 
Cn-l C n ,_2 C n _3 • • • Ci Cq 



(4.2) 
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Moreover, we set u n = exp(2ni/n) and put 



1 


1 


1 


1 


1 






, ,n— 1 


1 






2(n-l) 



1 w, 



ra-1 



2(n-l) 



Co' 



(n-l)(n-l) 



with % as the imaginary unit, and the matrix F n is called the Fourier matrix. Therefore, 
a circulant matrix can be diagonalized by the Fourier matrix F n , i.e., 



C n : = F n *diag(F n c)F n , 



(4.3) 



where diag(F n c) is a diagonal matrix holding the eigenvalues of C n . From (14. 3p . we can 
determine diag(F n c) in O(NlogN) operations by the FFT of the first column c of C n 

^.i. A O(NlogN) V-cycle MGM for ID 

We employ the V-cycle MGM to solve the one dimensional system (12. lip and illustrate 
the computational count of O(NlogN) per iteration and the required storage of O(N). 

Suppose A h = I — A k+1 / 2 , u h = U k+1 and f h = (I + A k+1 l 2 )U k + AtF fc+1 / 2 , then the 
resulting system (12. lip becomes the following general linear system 



A h u h = f h 



(4.4) 



and the system (I4.4p can be carry out by the Algorithm [T] ly, p. 443] and|2j 

In Algorithm [H at the highest (finest grid) level a mesh-size of h is used to solve the 
resulting system (I4.4p . The finest grid operator A^ is with the finest grid size h = Ax; 
the coarse grid operator An = 1 — A 2 ih, H = 2 l h, for 1 < I < log 2 iV — 1; h is the coarsest 
mesh-size; Ijj, Iff are respectively the prolongation operator and the restriction operator. 



For one dimensional system, the restriction operator/^ is defined by (16 1 

1 2 1 



rH 



1 

4 



1 2 1 



(4.5) 



2(lff) T . The smoothing operator smooth may be 



and the prolongation operator 1^ 
written as 

smooth^, m , f h ) = S h u + (7 - S^A^fh 



(4.6) 
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where Sh is the iteration matrix of the smoothing operator, and we define the weighted 
(damped) Jacobi iteration matrix by 0, p. 9] 



Shi 



I - coD~ 1 A h 



(4.7) 



with the weighting factor u G R, and D is the diagonal of Ah- Thus, the ( 14. 6ft becomes 
the following weighted Jacobi iteration 



smooth^, u , f h ) = S h ,ojU + tuD 1 f h . 



(4. 



In Algorithm [TJ the factors V\ and of smooth 1 ' 1 (Ah, uq, fh) and smooth^ 2 (Ah, Uh, fh) 
denote the number of weighted Jacobi iterations. In Algorithm [2J we give the stopping 
criterion of Algorithm [TJ 

Next, we illustrate the storage requirement of O(N) and the computational count of 
O(NlogN) per iteration. 

From (EUD, we have A k+1 / 2 = diag(£ fc+1 / 2 )A fc+1 / 2 , where 



^fc+1/2 



29? 9o+92 9s 
g« + g« 2g<* g« + g° g« 

g% + g2 2gf g% + g% 



9n x -2 



9n x -i 
9n x -2 



93 



9n x -2 
9n x -i 



9n x -2 



93 



9, 
9o+92 

2g? 



(4.9) 



being a Toeplitz matrix, and £ k+1 / 2 = [^ +1 ^ 2 , t^ 1 ^ ■ ■ ■ >£n 



2gf 
9S + 92 

k+l/2 ]T 
N x -ll ■ 

Then, we only need to store £ fc+1 / 2 anc [ g a = [2g±, g$ + g%, g$ . . . , g% x ^i\ T which have 
2N — 2 parameters, instead of the full matrix A k+1 l 2 which has (N — l) 2 parameters, i.e., 
the required storage O(N). Consider a one dimensional grid with N points, the finest 
grid, fl h , requires O(N) storage locations; Q 2h requires 2 _1 times as much storage as Q h ; 
Q 4h requires 4 _1 times as much storage as fl h ; in general, Q ph requires p~ l times as much 
storage as Q h . Adding these terms we obtain j^j 



Storage = O(N) 



1 1 1 

+ 2 + 22 + --'' + 2 lo s 2 ^-i 



0(N). 



Taking v a given vector, for the Toeplitz matrix vector multiplication A k+1 ^ 2 v, we first 
embed A k+l / 2 into a (2N X — 2) x (2N X — 2) circulant matrix, i.e., 



(4.10) 







V 




' A k+i/2 t ; 


* A k+l / 2 









t 
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Thus, using f!4.3j) . the computational count of A k+l / 2 v remains as O(NlogN). Therefore, 
for a V-cycle MGM, each level is visited O(NlogN) and grid VP h requires p work units. 
Similarly, adding these count we have 

V-cycle MGM computational count 

= O(NlogN) ■ ^1 + 1 + 1 + ..., +^ ZT ) = O(NlogN). 

By similar analysis, we know that the required storage is still O(N) and the compu- 
tational count O(NlogN) for multidimensional case. 

4.2. A O(NlogN) V-cycle LOD-MGM for 2D 
For D-AD scheme (Q-([33]), take 

A h x =I ~ A k+ D % u hx = U*; f hx = (J + 4+% + 2A k 2 + D %)U k + AtF™^ 
A hy = I-A k 2 + D % u hy =U^- h y = U*-A k 2 + D %U k . 

Similarly, for PR-AD scheme fl3.5l) - f)3.6p . denote 

A h x = I - A k +% u hx = U*; f hx = (J + A k Z> k + y F k + 1/2 ; 
A hy = I~A k+ D % u hy = U^- h y = {I + AtLl)U* + ^F^I\ 

Then, both the D-AD and PR-AD schemes, defined by (13.71) . reduce to the following LOD 
form: 

A hx u hx = f hx , (4.11) 

A h u hy = f hy . (4.12) 

Therefore, we can solve the two dimensional system (13 .7p by V-cycle LOD-MGM (see 
Appendix Algorithm [Q13]). 

The Algorithm [3] starts with the initial time t = and executes as follows: 

(1) First for every fixed y = yj (j = 1, . . . , N y — 1), using Algorithm [1] to solve a set of 

N x — 1 equations defined by (14. lip at the mesh points 2,;, i = 1, . . . , N x — 1, to get 

(2) Next alternating the spatial direction, and for each fixed x = Xi (i — 1, . . . , N x — 1) 

solving a set of — 1 equations defined by (I4.12p at the points yj, j = 1, . . . , N y — 1, 
once again we employ Algorithm [TJ to get uu ■ 
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4.3. A O(NlogN) V- cycle LOD MGM for 3D 
For D-AD scheme (I3.8I) - (I3.10I) . similarly, we set 

fh x = (/ + 4 + nll + + 2A k +f z )U k + AtF k ^ 2 - 

A, -r_A k+1/2 - 11, -n k ' 2 - f, - n** 1 - A k+1/2 n k - 

A -T A k+l l 2 - „ — TT k + 1 - f - n k > 2 A k+l / 2 TT k - 
A h z - 1 - A 3D,Az> U hz - U ) Jh z — U ~ A 3D,Az U ' 

Then, the D-AD scheme defined by (13. lip , becomes the following form: 

A hx u hx = f hx , (4.13) 

A hy u hy = f hy , (4.14) 
A hz u hz = f hz . (4.15) 
Therefore, we can solve the three dimensional system (13. lip by Algorithm [TJ |2]and|4l 

4-4- Convergence analysis 

In this subsection, we discuss the convergence of LOD-MGM. For the convenience 
of convergence analysis, we assume the coefficient c(x,y,z,t) is constant, and then all 
above are equal, and it can be denoted as £. Let's first consider the one 
dimensional case. 

From (14.91) . there exists 

A h = I- A k+l ' 2 = I - iA k+l ' 2 ee [cy-fcil^Dx^-i), (4.16) 

where 

c = 1 - 2g?t; d = -(g% + g^; c k = g% +1 £, k = 2, . . . , N x - 2, (4.17) 

and Ah is a symmetric Toeplitz matrix. From the proof of Theorem 3.2, we know that 
the matrix Ah is symmetric and strongly diagonally dominant with positive diagonally 
elements, then Ah is symmetric positive definite 0, p. 3]. 

Since the matrix Ah is symmetric positive define, we can define the following three 



different inner products [15j, p. 78] 



(u,v) = (Du,v), (u,v)i = (A h u,v), (u,v) 2 = (D A h u,A h v), (4.18) 

where D is the diagonal of Ah and along with their corresponding norms 1 1 • | |j {i = 0, 1, 2). 
If taking the coarsest grid size H = = 2h, Algorithm [1] is called the two-grid method 
(TGM). The TGM is rarely used in practice since the coarse grid operator may still be too 
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large to be solved exactly. However, it is useful from a theoretical point view as the first 

23| . Since 



12 



15 



16 



step to study the MGM convergence usually begins from the TGM 
;he MGM convergence analysis is still a challenge topic in computational mathematics 
3|. In the following we only consider the convergence of the TGM. 

Lemma 4.1. Let Ah be a M-matrix and the weighted factor < u < 1, then weighted 
Jacobi iteration \4-£ty converges. 

Proof. Taking M = D/u and N = D/oj — Ah, since < u < 1, then M and N be a 



regular splitting of a matrix Ah- Note that Ah is an M-matrix 
p. 119] 

p(S t J)=p{I-uD- x A h ) < 1. 



p. 3], thus we have 
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□ 



Lemma 4.2 



15 



p. 84]. Let Ah be a symmetric positive definite and rjo > p(D~ l Ah). 
If a < u(2 — oorjo), then the Jacobi relaxation with relaxation parameter < u < 2/t)q 
satisfies 



Sh,u,e h \\l < \\e h \\l - o-\\e h \\l, Ve^ 6 



>JV*-1 



(4.19) 



The inequality ( 14.191 ) is called the smoothing condition. For the TGM, the correction 



operator is given by [2|, p. 85] 



1 i 
E 



T c = I-I h H {A 



l T H A, 



•T c (Sh t 



h; see 
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, p. 89]. For 



therefore, the convergence factor of the TGM is 11(5^ 
convenience, we take v\ = and u 2 — 1. Therefore, the convergence factor of the TGM 
is given by H^^- T c ||i. 

Lemma 4.3. [15j, p. 89] Let Ah be a symmetric positive definite matrix and Sh,ui satisfy 



mm eh 

e H 6K iV x /2-l 



l[4.19\ ). Suppose that the interpolation l\ has full rank and that, for each eh, 

Ve h e R Nx ~\ (4.20) 
with k > independent of eh- Then, k > o and the convergence factor of the TGM 



jh 1 1 2 ^ II 1 1 2 



convergence factor satisfies \ \Sh,u • ^c||i < \J\ — o j k. 

Letting L^ x _\ = tridiag(— 1, 2, —1) be the (N x — 1) x (N x — 1) one dimensional discrete 
laplacian, then Ln x -i is a symmetric positive definite matrix. We define A rest = A h+1 ^ 2 + 
ciLjv^-i, where c\ is defined by (I4.17P and it can also be shown that A rest is symmetric 
and diagonally dominant with positive diagonally elements; then A rest is positive definite. 
Hence, we have the following equation 



(e h , A h e h ) = (e h , (I - c 1 L Nx ^ 1 + A rest )e h ) > {e h , (I - ciL,N x -i)eh), Ve^ G 



<>N X -1 



(4.2i; 
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Theorem 4.4. Since Ah, defined by \4-4ty j is a symmetric positive definite matrix, if 
taking a = u{2 — ur] ) with oj G (0,1], then Sh,u satisfies h4-19\ ) and the convergence 
factor of the TGM satisfies 



Sh,w ' T c \\i < y 1 — 2cr/5 < 1, (4.22) 



where rj = p(D~ l A h ) < 2. 



Proof. From Lemma 4.1, we have p(D l Ah) < 2. Taking r] and a in Lemma 4.2 as 
p{D~ x Ah) and a = u(2 — our) ), respectively, then Sh tU satisfies f)4.19p . 



Similar to the proof given in [l2j, we denote eh = (ei, e%, . . . , eN x -i)' T £ M^ -1 , eo = 
eN x = 0, and e# = (e 2 ,e 4 , . . . ,eN x -2) T £ R Nx//2 ~ l ; the norm || • || is defined by (14.181) . 
and D = diag(A^) = c I, where c is defined by (14.171) . There exists 



Nx/2-l, v 2 N x -l 

\\e h - IhZhWI = c ^2 ( e 2i+i - 2 e2i ~ 2 62i+2 ) ~ c °X^( e « 2 _ e * e *+0 ■ 

i=0 ^ ' t=l 

From the above inequality, we obtain 

iVs-l AT X — 1 

1=1 8=1 

Similarly, we can check that 

N x -1 N x -1 



$> 2 > - E^^ 1 - ( 4 - 23 ) 



t=i t=i 

Combining ( I4.2ip with ( 14.23(1 . we have 

N x -1 



lle^Hi = (e h ,A h e h ) > (e h , (I - ciL Nx _i)e h ) = ^((1 - 2ci)e 2 + 2c 1 e i e i+1 ) 

i=l 

Hence, Eq. (14301) holds, i.e., 

I |eft - ^r e i?||o < K l|e/v||i, 
c 2c 2 + 4£(-<) _ / 5\ 

since, according to Lemma 3.1, it is easy to check that g$ + g% < —g± < § (#0 + #2 ) ■ From 
Lemma 4.3, we have \\Sh,u • T c \\i < yl — 2cr/5 < 1. □ 
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Theorem 4.4 shows that the V-cycle scheme has a convergence factor with bound, that 
is independent of Ax. 

For the multidimensional case, we can similarly define the weighted Jacobi iteration 
matrix 

S hx , w = I -uD^A hx , S hyiW = I -uD^A hy , S hz , w = I -uD^Ah, 
and the LOD-TGM correction operators 

T Cx = I- I&Ah)- 1 !*^, T Cy = I- I&Ah)- 1 !*^, T Cz — I — i* (jI*)" 1 I» A hz , 



then the convergent factors of the LOD-TGM satisfy HS^^ • T Cx \\i < yj 1 — 2cr/5 < 1, 
H^-^Hi <y/l- 2a/5 < 1 and \\S hz>UJ ■ T Cz \\ x < yjl - 2a/5 < 1. 
Remark 4.7 As mentioned in the Introduction section, nowadays there are two differ- 
ent second-order discretization schemes for fractional operators 17|, |20j; all the analysis 



given in this paper can be parallel extended to the case that the fractional operators are 



discretized by the scheme given in 20J, in fact the scheme given in 20j has more wide 



applications because of its good p rop erties; and in Table HI we show the numerical results 
obtained by using the scheme of [20j to discretize the fractional operators of (12.131) . 

5. Numerical results 

We employ the V-cycle MGM and V-cycle LOD-MGM described in Section 4 to solve 
the one dimensional case (12.81) and multidimensional case (I2.13|ll.ip . respectively. The 
stopping criterion is taken as 

\\r il % r 

M 112 < 1Q - 7) 



|| r (°)||2 

where is the residual vector after I iterations. In all tables, N t denotes the number of 
time steps; N x , N y and N z , respectively, denotes the number of spatial grid points in x, 
y, and z direction, and the numerical errors are measured by the norm, 'Rate' denotes 
the convergent orders. 'CPU' denotes the total CPU time in seconds (s) or minutes (m) 
for solving the resulting discretized systems, and 'Iter' denotes the average number of 
iterations required to solve a general linear system A^Uh = fh at each time level. 

All numerical experiments are programmed in Python, and each computation was 
carried out on a PC with the configuration: AMD Phenom (tm) II X4 830 CPU 2.79 
GHZ and 3 GB RAM and a Linux operating system. All the numerical results listed 
in the following tables are got by the V-cycle MGM or V-cycle LOD-MGM with the 
parameters: the number of iterations (u 1 , v 2 ) = (1, 1) and (uj pre , uj p ost) = (1, 1/2). 
Remark 5.1. From our numerical experiences, we find that: 
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(1) With the increasing of the order of fractional derivative a from 1 to 2, the condition 

number of the matrix A k+1 ' 2 becomes bigger and bigger; and when getting the same 
accuracy the cost for a = 1.9 almost double the cost for a = 1.1; 

(2) For making MGM more efficient, the parameters can be dynamically chosen as when 

a increases from 1 to 2, correspondingly u post decreases from 1 to 0.5 and fixing 
u pre = 1 or fixing Upost = 1 and cu pre decreases from 1 to 0.5; 

(3) MGM is still powerful for second-order schemes, when simulating (12. 8p with the 

parameters given in subsection 5.1, in Table [T] it is shown that the second-order 
scheme costs 4.82 s to obtain the accuracy with the maximum error 1.2407e — 006, 



12[ costs 413.95 s when getting the accuracy with 



and the first-order scheme used in 
the maximum error 2.0358e — 006. 

5.1. Numerical results for ID 

Let us consider the one dimensional Riesz fractional diffusion equation (12.81) . where 
< x < 1 and < t < 1, with the variable coefficient c(x,t) = x a t, the forcing function 



f(x,t) = -e"V(l -x) 



+ 



x a te 



cos(aii/2) 



x 2 ~ a + (1 - xf~ a x 3 - a + (1 - x) 3 ~ a x 4 ~ a + (1 - x) 4 -° n 



- 6 —i A- + 12- 



r(3-a) I\4-a) r(5 - a) 

and the initial condition u(x,0) = x 2 (l — x) 2 and the boundary conditions u(0,t) = 
u(l,t) = 0. This fractional PDE has the exact value u(x,t) = e~*x 2 (l — x) 2 , which may 
be confirmed by applying the fractional differential equations 



r(p + i) 
r( P + i-v) 

r(p + i) 
r ^r - r(p + 1 _^ 



xl d v x {x - x L y = - ^ 1 A x - x L y-\ 



■MJxr - xf = - ^ ' A x R - xY-\ 



From Table [U we numerically confirm that the numerical scheme has second-order 
accuracy and the computational cost is of (^(iVlogiV) operations. 

5.2. Numerical results for 2D 

Consider the two dimensional Riesz fractional convection diffusion equation (I2.13p . on 
a finite domain 0<x<l,0<?/<l,0<t<l, and with the variable coefficients 

c(x, y, t) = x a y, d(x, y, t) = xy 13 , 

and the initial condition u(x, y, 0) = x 2 (l — x) 2 y 2 (l — y) 2 and the Dirichlet boundary con- 
ditions on the rectangle in the form u(0, y, t) = u(x, 0, t) — and u(l, y, t) = u(x, 1, t) — 
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Table 1: MGM to solve the resulting system (12. lip of the ID Riesz fractional convection 
diffusion equation (12. 8p at t — 1 and N t — N x . 



N t , N x 


a = 1.1 


Rate 


Iter 


CPU 


a = 1.9 


Rate 


Iter 


CPU 


2 5 


7.8755e-005 




4.0 


0.19 s 


7.5578e-005 




6.0 


0.28 s 


2 6 


2.1801e-005 


1.8530 


4.0 


0.45 s 


1.9255e-005 


1.9727 


6.0 


0.74 s 


2 7 


5.6999e-006 


1.9354 


4.0 


1.29 s 


4.8923e-006 


1.9766 


6.0 


2.14 s 


2 8 


1.4565e-006 


1.9684 


3.0 


2.56 s 


1.2407e-006 


1.9794 


6.0 


4.82 s 


2 9 


3.8540e-007 


1.9181 


3.0 


7.45 s 


3.1420e-007 


1.9814 


6.0 


13.56 s 


2 io 


9.7292e-008 


1.9860 


3.0 


18.63 s 


8.1028e-008 


1.9552 


6.0 


34.99 s 



for all t > 0. The exact solution to this two dimensional Riesz fractional convection 
diffusion equation is 

u(x, y, t) = e -*x 2 (l - x) V (1 - yf. 
From the above given quantities, it is easy to obtain the forcing function f(x,y,t). 



Table 2: LOD-MGM to solve the resulting system of the 2D Riesz fractional convection 
diffusion equation (12131) by the D-AD scheme (l2TT7j) - (l2TT8|) at t = 1 and N t = N x = N y . 



Nt, N x , N y 


a = 1.1, /3= 1.1 


Rate 


Iter 


CPU 


a = 1.8, p = 1.9 


Rate 


Iter 


CPU 


2 4 


2.4698e-005 




4.5 


2.09 s 


2.5475e-005 




7.0 


3.51 s 


2 5 


6.1249e-006 


2.0117 


4.0 


11.26 s 


6.5211e-006 


1.9659 


6.0 


18.03 s 


2 6 


1.5212e-006 


2.0095 


4.0 


55.85 s 


1.6662e-006 


1.9686 


6.0 


90.73 s 


2 7 


3.7812e-007 


2.0083 


4.0 


5 m 17 s 


4.2362e-007 


1.9757 


6.0 


8 m 45 s 


2 8 


9.4076e-008 


2.0070 


3.0 


21 m 36 s 


1.0744e-007 


1.9792 


6.0 


39 m 1 s 



From Table [2] and Table [3l numerically it can also be noticed that the D-AD and 
PR- AD are equivalent in two dimensional problems. We employ the LOD-MGM to solve 
two dimensional Riesz fractional diffusion equation, numerical results further display the 
computational cost is of 0(N\ogN) operations and the numerical scheme is second-order 
convergent. 

In particular, we further numerically confirm that this paper is still valid if the Riesz 
fractional derivative fll.2p is discretized by another existing second-order discretization 
scheme given in [20J], see Table HI In fact, theoretically we can also easily draw the same 
conclusion. 
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Table 3: LOD-MGM to solve the resulting system of the 2D Riesz fractional convection 
diffusion equation (1233)1 by the PR-AD scheme (1235]) - (j2~20ft at t = 1 and jV t = jV x = AT y . 



JV t) jV*, Ny 


a = 1.1, /3 = 1.1 


Rate 


Iter 


CPU 


a = 1.8, P = 1.9 


Rate 


Iter 


CPU 


2 4 


2.4698e-005 




4.5 


2.05 s 


2.5475e-005 




7.0 


3.48 s 


2 5 


6.1249e-006 


2.0117 


4.0 


11.17 s 


6.5211e-006 


1.9659 


6.0 


17.71 s 


2 6 


1.5212e-006 


2.0095 


4.0 


55.14 s 


1.6662e-006 


1.9686 


6.0 


90.24 s 


2 7 


3.7812e-007 


2.0083 


4.0 


5 m 18 s 


4.2362e-007 


1.9757 


6.0 


8 m 42 s 


2 8 


9.4075e-008 


2.0070 


4.0 


22 m 25 s 


1.0744e-007 


1.9792 


6.0 


39 m 31 s 



Table 4: LOD-MGM to solve the resulting system of the 2D Riesz fractional convection 
diffusion equation fl2TT3|) by the D-AD scheme ll2~T71) - (l2~T8l) at t = 1 and N t = N x = PL; 
the Riesz fractional derivative (11.21) is discretized by the scheme given in [20|]. 



N t , N x , Ny 


a = 1.1, p = 1.1 


Rate 


Iter 


CPU 


a = 1.8, P = 1.9 


Rate 


Iter 


CPU 


2 4 


2.4592e-005 




5.0 


2.25 s 


2.4532e-005 




7.0 


3.5 s 


2 5 


6.1745e-006 


1.9938 


4.0 


11.42 s 


6.0897e-006 


2.0102 


6.0 


17.98 s 


2 6 


1.5426e-006 


2.0010 


4.0 


56.57 s 


1.5102e-006 


2.0116 


6.0 


90.74 s 


2 7 


3.8444e-007 


2.0045 


4.0 


5 m 28 s 


3.7350e-007 


2.0155 


6.0 


8 m 41 s 


2 8 


9.5743e-008 


2.0055 


4.0 


22 m 43 s 


9.2374e-008 


2.0155 


6.0 


38 m 22 s 



5. 3. Numerical results for 3D 

Consider the three dimensional Riesz fractional convection diffusion equation (jl.ip . on 
a finite domain < x < 1, < y < 1, < z < 1, and < t < 1, and with the variable 
coefficients 

c(x, y, z, t) = x a yz, d(x,y,z,t) = xjfz, e(x,y,z,t) = xyz r , 

and the initial condition u(x, y, z, 0) = x 2 (l — x) 2 y 2 (l — y) 2 z 2 (l — z) 2 and the zero Dirichlet 
boundary conditions on the cube. The exact solution to this three dimensional Riesz 
fractional convection diffusion equation is 

u(x, y, z, t) = e-V(l - xfy 2 {\ - y) 2 z 2 (l - z) 2 . 

According to the above conditions, it is easy to obtain the forcing function f(x,y,z,t). 

The numerical results in Table [5] are obtained by employing the LOD-MGM to solve 
the three dimensional Riesz fractional diffusion equation, they again display the compu- 
tational count of O(NlogN) operations and the scheme is second-order convergent. 
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Table 5: LOD-MGM to solve the scheme fl2T23|) -f l2T25j) of the 3D Riesz fractional convec- 
tion diffusion equation ( II. ip at t — 1 and N t = N x = N y = N z . 



N t 


a = P = 7= 1.1 


Rate 


Iter 


CPU 


a = 1.8, P= 1.9,7 = 1-8 


Rate 


Iter 


CPU 


2 3 


5.9349e-006 




4.75 


3.31 s 


5.8311e-006 




7.0 


5.28 s 


2 4 


1.4792e-006 


2.0044 


4.0 


42.19 s 


1.4867e-006 


1.9717 


7.0 


70.1 s 


2 5 


3.7377e-007 


1.9846 


4.5 


7 m 40 s 


3.8428e-007 


1.9519 


6.0 


13 m 1 s 


2 6 


9.3376e-008 


2.0010 


4.37 


76 m 13 s 


9.8179e-008 


1.9687 


6.0 


136 mis 



6. Conclusions 

With the appearing of the two effective second-order discretization schemes [l?], 0] 
and MGM being successfully employ to solve the resulting system of the one dimensional 
fractional diffusion equation discretized by first-order scheme 12[ , our attentions turn to 
the possibility of efficiently solving the multidimensional fractional Riesz diffusion equa- 
tion by second-order scheme and MGM. This paper shows that when solving (II. ip the 
second-order schemes are unconditionally stable, and the structure of the resulting ma- 
trix algebraic equations is almost the same as the one by the first-order scheme and then 
none computational costs increase but accuracy is greatly improved if using the second 
scheme instead of the first-order one. And LOD-MGM still preserves its powerfulness of 
0(N\ogN) computational counts and of O(N) storage when solving the resulting ma- 
trix system of the multidimensional fractional Riesz diffusion equation discretized by the 
second-order scheme. For computing the completely same equation, from Table HJ it can 
be noticed that the second-order scheme costs 4.82 s to obtain the accuracy with the 
maximum error 1.2407e — 006, and the first-order scheme used in 12[ costs 413.95 s when 
getting the accuracy with the maximum error 2.0358e — 006. 

Last but not least, we want to refer to that although this paper focus on using the 
second-order discretization given in 17J, all the analysis of this paper is still valid if 
applying the second-order discretization in [20] ; in fact the validness is already verified 
numerically in Table HI 
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Appendix 

For a general linear system 

A h u h = f h , 

we employ the following V-cycle MGM (Algorithm |TH2j) to solve one dimensional ( 12. 81) 
and V-cycle LOD-MGM (Algorithm [TH) solve two dimensional (I2.13p . Solve the three- 
dimensional system (11. ip by Algorithm [Tf2l and HJ 

Algorithm 1 MGM for ID u h =V-cyc\e(A h , u , f h ) 
1: Pre-smooth: Uh := smooth 1 ' 1 (Ah, uq, fh) 
2: Get residual: r h = f h ~ A h u h 
3: Coarsen: th = Ih r h 
4: if H == ho then 
5: Solve: Ah^h = rn 
6: else 

7: Recursion: £ H = V-cycle(A^, 0, r H ) 
8: end if 

9: Correct: u h -=u h + Ijj^H 
10: Post-smooth: Uh := smooth." 2 (Ah, Uh, fh) 
11: Return Uh 



Algorithm 2 Stopping criterion for MGM 
1: u h := u 

2: r0 := \\A h u - f h \\2 

3: r h := r 

4: while — > e do 

ro 

5: u h := V-cycle(v4 h , u h , f H ) 
6: r h := \ \A h u h - f h \\ 2 
7: end while 
8: Return Uh 
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